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Abstract — When recovering a sparse signal from noisy com- 
pressive linear measurements, the distribution of the signal's non- 
zero coefficients can have a profound affect on recovery mean- 
squared error (MSE). If this distribution was apriori known, 
then one could use computationally efficient approximate message 
passing (AMP) techniques for nearly minimum MSE (MMSE) 
recovery. In practice, though, the distribution is unknown, moti- 
vating the use of robust algorithms like Lasso — which is nearly 
minimax optimal — at the cost of significantly larger MSE for 
non-least-favorable distributions. As an alternative, we propose 
an empirical-Bayesian technique that simultaneously learns the 
signal distribution while MMSE-recovering the signal — according 
to the learned distribution — using AMP. In particular, we model 
the non-zero distribution as a Gaussian mixture, and learn its 
parameters through expectation maximization, using AMP to 
implement the expectation step. Numerical experiments on a wide 
range of signal classes confirm the state-of-the-art performance 
of our approach, in both reconstruction error and runtime, in 
the high-dimensional regime. 



I. Introduction 

We consider estimating a A'-sparse (or compressible) sig- 
nal x <G M. N from M < N linear measurements y = 
Ax + w € R M , where A is known and w is additive 
white Gaussian noise (AWGN). For this problem, accurate 
(relative to the noise variance) signal recovery is known to 
be possible with polynomial-complexity algorithms when x 
is sufficiently sparse and when A satisfies certain restricted 
isometry properties 0, or when A is large with i.i.d random 
entries as discussed below. 

Lasso (or, equivalently, Basis Pursuit Denoising Q), is 
a well-known approach to the sparse-signal recovery problem 
that solves the convex problem 



Alasso = argmin \\y - Ax\\^ + A| asso ||±||i, 



(1) 



with A| asso a tuning parameter that trades between the spar- 
sity and measurement-fidelity of the solution. When A is 
constructed from i.i.d entries, the performance of Lasso can 
be sharply characterized in the large system limit (i.e., as 
K, M, N — s- oo with fixed undersampling ratio M/N and 
sparsity ratio K/M) using the so-called phase transition curve 
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(PTC) 0, (8). When the observations are noiseless, the PTC 
bisects the M/N-vstsus-K/M plane into the region where 
Lasso reconstructs the signal perfectly (with high probability) 
and the region where it does not. (See Figs. |2]|4]) When the 
observations are noisy, the same PTC bisects the plane into 
the regions where Lasso's noise sensitivity (i.e., the ratio of 
estimation-error power to measurement-noise power under the 
worst-case signal distribution) is either finite or infinite |9). 
An important fact about Lasso's noiseless PTC is that it is 
invariant to the distribution of the nonzero signal coefficients. 
In other words, if the vector x is drawn i.i.d from the pdf 



Px(x)=\fx(x) + (l-X)S(x), 



(2) 



where S(-) is the Dirac delta, fx(') is the active-coefficient 
pdf (with zero probability mass at x = 0), and A = K/N, 
then the Lasso PTC is invariant to fx{ )- While this implies 
that Lasso is robust to "difficult" instances of fx{ ), it also 
implies that Lasso cannot benefit from the case that /x(-) is 
an "easy" distribution. For example, when the signal is known 
apriori to be nonnegative, polynomial-complexity algorithms 
exist with PTCs that are better than Lasso's iflOll . 

At the other end of the spectrum is minimum mean-squared 
error (MMSE)-optimal signal recovery under known marginal 
pdfs of the form (|2|i and known noise variance. The PTC 
of MMSE recovery has been recently characterized fiTl and 
shown to be well above that of Lasso. In particular, for any 
fx(-), the PTC on the M /N-versus-K/M plane reduces to the 
line K/M = 1 in both the noiseless and noisy cases. More- 
over, efficient algorithms for approximate MMSE-recovery 
have been proposed, such as the Bayesian version of Donoho, 
Maleki, and Montanari's approximate message passing (AMP) 
algorithm from Ifl2l . which performs loopy belief -propagation 
on the underlying factor graph using central-limit-theorem 
approximations that become exact in the large-system limit 
under i.i.d A. Although AMP's complexity is remarkably 
low (e.g., dominated by one application of A and A J per 
iteration with typically < 50 iterations to convergence), it 
offers rigorous performance guarantees in the large-system 
limit [13]. To handle arbitrary noise distributions and a wider 
class of matrices A, Rangan proposed a generalized AMP 
(GAMP) lfl4l that forms the starting point of this work. (See 
Table H) 

In practice, one ideally wants a recovery algorithm that does 
not need to know px(-) and the noise variance a priori, yet 
offers performance on par with MMSE recovery, which (by 
definition) requires knowing these prior statistics. Towards this 
goal, we propose a recovery scheme that aims to learn the prior 
signal distribution px(-), as well as the variance of the AWGN, 
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while simultaneously recovering the signal vector x from the 
noisy compressed measurements y. To do so, we model the 
active component fx(-) in © using a generic L-term Gaussian 
mixture (GM) and then learn the GM parameters and noise 
variance using the expectation-maximization (EM) algorithm 
[15]. As we will see, all of the quantities needed for the 
EM updates are already computed by the GAMP algorithm, 
making the overall process very computationally efficient. 
Moreover, GAMP provides approximately MMSE estimates of 
x that suffice for signal recovery, as well as posterior activity 
probabilities that suffice for support recovery. 

Since, in our approach, the prior pdf parameters are treated 
as deterministic unknowns, our proposed EM-GM-AMP algo- 
rithm can be classified as an "empirical-Bayesian" approach 
lfl6l . Compared with previously proposed empirical-Bayesian 
approaches to compressive sensing (e.g., lfT7l - lfT9l ), ours has 
a more flexible signal model, and thus is able to better match a 
wide range of signal pdfs px{'), as we demonstrate through a 
detailed numerical study. In addition, the complexity scaling of 
our algorithm is superior to that in lfT71 - lfT9l , implying lower 
complexity in the high dimensional regime, as we confirm 
numerically. 

Notation: For matrices, we use boldface capital letters like 
A, and we use tr(A) and to denote the trace and 

Frobenius norm, respectively. Moreover, we use (-) T to denote 
transpose, (•)* conjugate, and (-) H conjugate transpose. For 
vectors, we use boldface small letters like x, and we use 
\\x\\ p = (^2 n |x„| p ) 1//p to denote the £ p norm, with x„ 
representing the n th element of x. For a Gaussian random 
vector x with mean m and covariance matrix Q, we denote 
the pdf by M(x; m, Q), and for its circular complex Gaussian 
counterpart, we use CAf(x;m,Q). Finally, we denote the 
expectation operation by E{-}, the Dirac delta by 5(-), the 
real field by R, and the complex field by C. 

II. Gaussian-Mixture GAMP 

We first introduce Gaussian-mixture (GM) GAMP, a key 
component of our overall approach, where the coefficients in 
x = [xi, . . . , xjv] t are assumed to be i.i.d with marginal pdf 



L 

p x {x;X,oj,e,4>) = (l-A)5(x)+A^o;^(x;^,^), (3) 

=i 



where S(-) is the Dirac delta, A is the sparsity rate, and, 
for the k th GM component, ujf., 9k, and fa are the weight, 
mean, and variance, respectively. In the sequel, we use u> = 
[oji, . . . ,ojl] J and similar definitions for and <f>. The noise 
w = [wi, . . . , wm ] T is then assumed to be i.i.d Gaussian, with 
mean zero and variance tp, i.e., 



p w (w;ip) = Af(w;0,ip), 



(4) 



and independent of x. Although above and in the sequel we 
assume real-valued quantities, all expressions in the sequel can 
be converted to the circular-complex case by replacing Af with 
CAf and removing all i's. We note that, from the perspective 
of GM-GAMP, the prior parameters q = [A, u), d, <f>, tp] and 
the number of mixture components, L, are treated as fixed and 
known. 



GAMP models the relationship between the m th observed 
output y m and the corresponding noiseless output z m = al n x, 
where aJ m denotes the m th row of A, using an arbitrar 
conditional pdf PY\z(ym\z m )- Under the AWGN assumptior 1 
©' Py\z(u\ z ) — N{y; z, tp), and thus the quantities g u\ 
and g ou t( ) in Table [T] become lfT4l 



9out{y,z,fi z ;q) = 



y- z 



, ffout(y^^ z ;g) = 



(5) 



As for the quantities g m {-) and g' m {-) in TableU straightforward 
calculations with our GM signal prior (0) reveal that 



g- m (r,fi r ;q) = Tr(r, t i r ;q) 

V r g' n (r,H r \q) = -\g\n{?,n\q)\ 2 + 7r(f , <?) 

J2e=i Pt{r, /i r ; q) (\li{r, /i r ; q)\ 2 + v t {f, //'; q)) 



where 



J2i=iMr : n r ;q) 



{ ** r \ A 

n^A* ;g) = - 



jt{f,fi r ;q) 



a ?/n r + 9 e /<f>t 



v i\r , A* 5 q) = 



l/fl r + l/4>£ 

1 



(6) 

(7) 

(8) 
(9) 

(10) 
(11) 



1/^ + 1/0/ 

TableU then implies that GM-GAMP's marginal posteriors are 

Px(x n ;q) N{x n ;f n ,n r n ) 



p(x n \y;q) = 



((r n ,n r n ;q) 



1 



C(f,A* r ;9)= px(x;q)M(x;f,n r ) 



(12) 

Af(x n ;r n ,n r n ) 

CC^n.A*^;?) 

(13) 
(14) 



= (l-A)AA(O;f, / i r ) + A^^AA(O;f-c9 £ , / i r + 0,). (15) 



From (fT3T l. it is straightforward to show that the posterior 
support probabilities returned by GM-GAMP are 



Pr{x„ ^ | y; q} = n(f n ,fj^; q) 



(16) 



In principle, one could specify GAMP for an arbitrary signal 
prior px(-)- However, if the integrals in (D4)-(D6) are not 
computable in closed form (e.g., when px(') is Student's-t), 
then they would need to be computed numerically, thereby 
drastically increasing the computational complexity of GAMP. 
In contrast, for GM signal models, we see from the above that 
all steps can be computed in closed form. Thus, a practical 

'Because GAMP can handle an arbitrary py\z ('!')> me extension of EM- 
GM-AMP to additive non-Gaussian noise, and even non-additive measurement 
channels (such as with quantized outputs or logistic regression), is straightfor- 
ward. Moreover, the parameters of the pdf p Y | z ( ' I ' ) could be learned using a 
method similar to that which we propose for learning the AWGN variance -0 
in Section Hn] Finally, one could even model Py|z(|) as a Gaussian mixture 
and learn the corresponding parameters. 
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definitions: 

Vz\y{Av\~ z ^ z ) 

Sin(r,/0 

initialize: 

Vn : x n (l) 
Vn : n£(l) 
Vm : u m (0) 
for t = 1 : T max , 
Vm : z m (t) 
Vm : /^(t) 
Vm : pm (t) 
Vm : « m (f.) 
Vm : /CM 
Vn :«£(*) 
Vn : r„ (t) 
Vn:u£(t+1) 
Vn : x»(t+l) 



PY Af(z;z,j^) 

Pi-|z(al 2 ') Af(z';i,M*) 

^(EjriH'ltfi*,/*"}-*) 
1 / var z|y {z|y;z,^"} 

px(^)-tf(j;' ! .<' r ) 
Px( l ')Af( l '; f ,i''') 

Jp*Px|Y<*lw;<'iAt r ) 



77=- /a I 1 - Sin(r, M r )l 2 Px|y(^|y; r, « r ) 



4 |x-£ n (l)| 2 p x (i) 


E (t) 

E„=l klm«| 2 ^(*) 
«m(*) - «m(*) «m(* ~ 1) 

(*),*4,(*)) 
-9' 0Ui {ym,p m {t), [l Z m {t)) 

(ELl |Am»| a ^(*)) _1 
Ei=l^„ 

„(*»(*)>/£(*)) 

Sin (*)>/"£.(*)) 



.(*) 



end 



ifE»=l |*n(*+l) - £n(t)| 2 < r gamp ELl IM*)I 



(Dl) 
(D2) 
(D3) 
(D4) 
(D5) 
(D6) 

(ID 
(12) 
(13) 



(Rl) 
(R2) 
(R3) 
(R4) 
(R5) 
(R6) 
(R7) 
(R8) 
(R9) 

, break (RIO) 



TABLE I 
The GAMP Algorithm (14) 

approach to the use of GAMP with an intractable signal prior 
Px(-) is to approximate px{ ) using an L-term GM, after 
which all GAMP steps can be easily implemented. The same 
approach could also be used to ease the implementation of 
intractable noise priors py\z (■]■)■ 

III. EM Learning of the Prior Parameters q 
In our approach, the expectation-maximization (EM) algo- 
rithm |[T5l is employed to learn the prior parameters q = 
[X,u>,0,<p,ip], The EM algorithm is an iterative technique 
that increases a lower bound on the likelihood p(y, q) at each 
iteration, thus guaranteeing that the likelihood converges to a 
local maximum. In our case, the "hidden data" is chosen as 
{x, w}, implying the iteration-i EM update 

q l+1 = argmaxE { \np(x, w; q) I y; q 1 }, (17) 

where E{-\y;q 1 } denotes expectation conditioned on the ob- 
servations y under the parameter hypothesis q % . Since it is 
impractical to update the entire vector q at once, we update 
q one component at a time (while holding the others fixed), 
which can be recognized as the "incremental" technique from 
J20). In the sequel, we use "Q\ A " to denote the vector q % with 
the element A removed (and similar for the other parameters). 
We emphasize that our EM-guided GM-fitting procedure dif- 
fers from the standard one (e.g., ||2T1 p. 435]) because, in our 
case, realizations of x are not directly observed. 

A. EM Update of the Gaussian Noise Variance tp 

We first derive the EM update for the noise variance 
if) given a previous parameter estimate q % . Because w is 
apriori independent of x and i.i.d, we can write p(x, w; q) = 
CY\m = iPw{wm'i4>) f° r a "0 -invariant constant C, and so 

M 

ip t+1 =argmaxV E { \np w (w m ; ip) \ y 1 q 1 }. (18) 



The maximizing value of ijj in (fT8l is necessarily a value of 
ijj that zeroes the derivative of the sum, i.e., that satisfies 

M , 

E / p( w m\y;q t )-Tj'^i-pw(w m ;ip) = o. (19) 

m=l )w ™ V 

Because pw(w m ', 4>) = ■A/"(Wm> 0, ip), it is readily seen that 



-lnp w ( Wm ^) = -l — -- J , 



(20) 



which, when plugged into ( fT9] l, yields the unique solution 

M 

r +1 = TjY] / \w m \ 2 p(w m \y;q l ). (21) 



Since w m = y m — z m for z m = a^x, we can also write 
1 M f 

mJ2 I \Vm- z m \ 2 p{z m \y;q l ) (22) 



V-' 



i+i 



m— 1 
M 



(23) 



where z m and p, z m , the posterior mean and variance of z m , are 
available from GAMP (see steps (R1)-(R2) in Table Hi. 

B. EM Updates of the Signal Parameters: BG Case 

Suppose that the signal distribution px(') is modeled using 
an L = 1-term GM, i.e., a Bernoulli-Gaussian (BG) pdf@ In 
this case, the marginal signal prior in © reduces to 



p x (x; X, uj, 6, 4>) = (1 - X)S(x) + XAf(x; 0, < 



(24) 



Note that, in the BG case, the mixture weight u is, by 
definition, unity and does not need to be learned. 

We now derive the EM update for A given previous 
parameters q l = [A', 9 l 1 (f> 1 , ip 1 ]. Because x is apriori in- 
dependent of w and i.i.d, we can write p(x,w;q) = 
CYln=iP x ( Xn 'i ^' ^' ^) ^ or a A-invariant constant C, and so 

N 

= argma x yE{lnpx(in;A,gu) \y,q 1 }- (25) 

Ae(o,D „ = i 

The maximizing value of A in ( 1251 ) is necessarily a value of A 
that zeroes the derivative of the sum, i.e., that satisfies 



N 

E 

71=1 



f d 
J p(xn\y,q l )-^lnp x (x n ;X,q{ x ) ^0. (26) 



For the BG px{x n ] A, 0, 

d , i , JV(. 
lnpx(x„; A, q XA ) = 



in (l24l . it is readily seen that 

9\ </>')- S(x n ) J 4- x n ^0 



Px(x„; A, q\ x ) 



X n = 0. 



(27) 



Plugging d27l i and ( fT3l into ( |26] |. it becomes evident that 
the neighborhood around the point x n = should be treated 
differently than the remainder of K. Thus, we define the closed 

2 We note that, after presenting our initial work on EM-BG- AMP in [T] 
and submitting it to 0, a closely related approach appeared in 1221 in the 
context of "seeded belief propagation." Although the main contribution of our 
current work is EM-GM-AMP, we detail the case of EM-BG-AMP in order to 
demonstrate that the approximations needed for the EM-GM recursion become 
tight in the L = 1 case (i.e., EM-BG). 
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ball B e = [— e, e] and its complement B e = K \ B e , and note 
that, in the limit £4 0, the following is equivalent to (l26l i: 

1 N I N f 

jJ2 _p{x n \y;g l ) = j—^ / ■ 



— 1 3!n £ 



l-7r(f„,/4;<f) 



(28) 



To verify that the left integral converges to the 7r(f n ,/i^; <jr l ) 
defined in ©, it suffices to plug : (TT~3T> into (l28T > and apply the 
Gaussian-pdf multiplication rulejj The right integral must then 
equal one minus the value of the left integral because, together, 
their regions of integration cover the real line. Finally, the EM 
update for A is the unique value satisfying (l28l l as e — > 0, 
which is readily shown to be 



1 N 

A i+1 = ^5>(r„,/4;<f 



(29) 



Conveniently, the quantities {n(r n , fi^; q % )}^—i are GM- 
GAMP outputs, as per ( [T6l l and ([9). 

Similar to ( |25T l, the EM update for 9 can be written as 



JV 



f +1 = argmax^E{lnp x (x„;0,g( )|y;< Z i }. (30) 



we 



Tl=l 



The maximizing value of 6> in (|30l l is again a necessarily a 
value of 9 that zeroes the derivative, i.e., that satisfies 

N f d 

J2j P(x n \y;q l )-^\np x (x n ;9,q\ e ) = 0. (31) 



For the BG px {x n ; A, 6, 4>) given in ( l24b . 

(x n -0) JW(i„;8,f) 



x„ # 

£„ = 0.' 



(321 



Splitting the domain of integration in ( f3TT > into £> e and S £ as 
before, and then plugging in (l32l i. we find that the following 
is equivalent to PTT ) in the limit of 6 ^ 0: 



V/ _(x n -9)p(x n \y:q l ) = 0. 
The unique value of 6> satisfying d33l as e 4 is then 



Qi+l 



1 N 

= vi + i W E ,r ( f "'''" ; ^ i (' : "'<'" ;g ') 



(33) 



(34) 



(35) 



for the GM-GAMP outputs {71 (f„, q 4 )}«=i defined in 
(fTOb . The equality in ( f33T > can be verified by plugging the 
GAMP posterior expression from ([T3l into d34T l and simplify- 
ing via the Gaussian-pdf multiplication rule. 

Similar to (IZ5l l. the EM update for <fr can be written as 



JV 



1 = argmax^E{lnj3x(a; n ;0,Q^) | y;q 1 }. (36) 



The maximizing value of <f> in ([36l l is again necessarily a value 
of <fi that zeroes the derivative, i.e., that satisfies 

N r d 

^2 P{xn\y]q l )-r-\-npx{x n ]4),q\^) ^Q. (37) 

n=l Jx « 9 

For the px(x n ', A, 9, 4>) given in (l24l . it is readily seen that 

(/ 



lnpx(a; n ; A\< 



1 ( \x n -8*\ 2 _ l\ A' Af(x n ; 9 i , <f>) 
2\ (0)2 cj>) p x {x n ;,4>,q\^ 

1 f _ 1 \ 2-^0 



X n = 



(38) 



Splitting the domain of integration in ( 137] ) into S e and B t as 
before, and then plugging in ( 1381 ), we find that the following 
is equivalent to ( f37T > in the limit of e ^ 0: 



JV 



1 



The unique value of <f) satisfying d39l as e ^ is then 

i+1 _ T,n=i l ™^o J Xn<£ B- t \x n - 9 i \ 2 p{x n \y;q l ) 



ip(x„|y;g l ) = 0. (39) 



(40) 



Finally, we expand \x n - 9 l \ 2 = \x n \ 2 - 2Re(x* n 6 z ) + l^l 2 
which gives 

1 N 

<P l+1 = X i +1N w (^n>M^>g*)(|0*-7i(rn,M^;a*)| 2 + Mrn,/4;a*)) 

(41) 

for the GM-GAMP outputs {vi(r n , if n ; 9*)}^Li defined in 
(fTTl i. The equality in fiTI ) can be verified by plugging the 
GAMP posterior expression from ([TBI into (|40l and simplify- 
ing using the Gaussian-pdf multiplication rule. 



C. EM Updates of the Signal Parameters: GM Case 

We now generalize the EM updates derived in Section ITlI-BI 
to the GM prior given in (f3]i for L > 1. As we shall see, it 
is not possible to write the exact EM updates in closed-form 

when L > 1, and so some approximations will be made. 

We begin by deriving the EM update for A given the 
previous parameters q l = [A\ u>\ 8 1 , (t> , jJj 1 ]. The first two 
steps are identical to the steps ( I25l > and ( 1261 presented for the 
BG case, and for brevity we do not repeat them here. In the 
third step, use of the GM prior (01 yields 



d , , . i . Eti"iAf(i„;4^-i(x„) 
~ lllPX(X " ; ^ = Px(^;A, qU ) ^ 



iJX 



Af{x;a,A)Af(x;b,B)=N{x; 1/aXi/b ' i/A+i/a W(0; a-b, A+B). update for A is that given in (^J 



x„ + 
x n = 
(42) 

which coincides with the corresponding BG expression ( l27l ). 
The remaining steps also coincide with those in the BG case, 
and so the final EM update for A, in the case of a GM0 is 
given by ( |29l . 

We next derive the EM updates for the Gaussian-Mixture 
parameters u>, 6, and <fc. For each k = 1, . . . , L, we incremen- 
tally update 9k, then (f>k, and then the entire vector u>, while 

4 The arguments in this section reveal that, under signal priors of the form 
px{x) = (1 — X)S(x) + Xfx(x), where fx{') can be arbitrary, the EM 
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holding all other parameters fixed. The EM updates are thus 



JY 



6\ +1 = argmaxVE{hip x (.T„;^,qV)ly;^}' < 43 ) 

N 

= argmax^EjlnpA^^^fe^^Jlyjq 4 } (44) 



71=1 



N 



u> l+1 = argmax ^ E { \np x (x n ; to, q\ u )\y, q*}{45) 

ui>0:$~J fc w fc =l n=1 

Following (13 It . the maximizing value of 6k in d43l is again 
necessarily a value of that zeros the derivative, i.e., 

N r d 

Xjy P(^n|y;9 i )^lnp Jf (a;„;0fc,q* Ne J = 0, (46) 

where p(x n \y;q l ) = px(x n ]q i )Af(xp:f n ,fjr n )/C(f„,fJ, r n -,q i ) 
from (D4), with GM px(x;q) from ([3]) and ((r,{i r ;q) from 
< Ti~5T >. Taking the derivative, we find 



we find that better reconstruction performance is obtained by 
fixing the means at zero (i.e., 6\ = Vft, i). Thus, in the 
remainder of the paper, we consider two modes of operation: a 
"sparse" mode where 6 is learned via the above EM procedure, 
and a "heavy-tailed" mode that fixes 6 = 0. 

Following d46l l. the maximizing value of <fik in d44l is 
necessarily a value of 4>k that zeroes the derivative, i.e., 

N f d 

E J p( x n\y\ql-^^px(x n -,(l)k,q\^ k ) = o. (52) 



(53) 



,10,, 



lnpx(,x n ;0 k ,q\B k ) 



As for the derivative in the previous expression, we find 

d t i f \x n -ej\ 2 i 

(1 - \*)5{x n ) + X'(ivlAf(x n ; 0%, cf> k ) + Y.^ k u\M(x n \ 0\, 4>\)) 

Integrating ( l52l separately over B e and B e , as in d28l i. and 
taking e —> 0, we find that the B t portion vanishes, giving 



(1 - \')S(x n ) + A'(w*A^(a;„;e fc) 1 fc ) + T, e , k ^Af(x n ; 6' e , ^ e )) 



(47) 



y f ES: 



p(i»|»n # 0, y; q , )A , tjJ.AA(x„; 9' k ,rp k ) 



)(a»«JV(x n ; 0j, + "iJV(a;„; 8j, 0j)) V 0* 



(54) 



Integrating d46b separately over S e and S e , as in d28l ), and 
taking e — > 0, we find that the S e portion vanishes, giving the 
necessary condition 

p(x n \x n ^ 0, y;q l )X l u} l k M{x n ;6 k ,(j> k ){x n - 9 k ) 



Similar to d48l i. this integral is difficult to evaluate, and 
so we again apply the approximation Af(x n ;6 k ,(j)k) ~ 
Af(x n ; 6 k , 4> l k ) in the numerator and denominator, after which 
several terms cancel, yielding the necessary condition 



^ r p(x n \x n ^ 0,y;q\ 

^[ Jxn arn,^ n ,q l )HM(x n ; 



N 



AT(o; n ;r 7i ,^)A i a;» ; Ar(a; n ;^, 9 ij ; ) ( \x n - 6 



(48) 

Since this integral cannot be evaluated in closed form, we 
apply the approximation J\[(x n ;6k,4> k ) ~ N{x n ,d l k ,(j) l k ) 
in both the numerator and denominator, and subse- 
quently exploit the fact that p(x n \x n ^ 0,y;q l ) = 

■A' (x n j T'n j 

Mn) He w\N(x n ;6\, (f>\) to cancel terms, giving the 
(approximated) necessary condition 

V- f ywjN (x n ;f n , [i r n )M (x n ; 6\ , 0| ) 

^1 0^7) (X " " " fc) = °- 

(49) 

We then simplify d49l using the Gaussian-pdf multiplication 
rule, and set k +1 equal to the value of 6k that satisfies d49l >: 



•Pk 



I =0. 

(55) 



To find the value of (f> k satisfying ( |55l l, we expand \x 



Qi 12 



2 Re« 61) + \6 k \ 2 and apply the Gaussian-pdf 
multiplication rule, which gives 

+ l_En = lPr{^0,^=fc|«;q i }(|el -7fc(^,^;g i )| 2 +»/ fc (f„,^;q)) 



(56) 



for Pr{a;„ 7^ 0, fc„ = fc|y; q 1 } defined in ( Bil l. 

Finally, the value of the positive uj maximizing d45l > under 
the pmf constraint J2k=i w fc = 1 can t> e found by solving the 
unconstrained optimization problem max u ^ J(<*>, £)> wnere £ 
is a Lagrange multiplier and 



JV 



ai+l _ E„=i Pr { ;r « 7^ 0, fc n = ft I y; q J }7fc(f„, q l ) 
1 k — 



En=i P Hx n ^ 0, k„ = k | y; q 4 } 



where we define 



(50) 



Pr{x n 7^ 0, ftn = k | y; q 1 } = 7r(f„, /j^, q* 



/3 fc (f tl ,/i;;q') 



71=1 ^£=1 

(57) 

^ /" ( L 

= Y1 / K^li;;9')inpx(j:n;w,qy-( Xl Wf_1 

n=l ^ £=1 



(51) 

using /3fc(f, /i r ; q 8 ) from (O and 7fe(f, q l ) from ( fTOb . Here, 
the notation "fc n = ft" describes the event that x n was generated 
from mixture component ft. 

For sparse signals x, we find that learning the GM means 
{6k} using the above EM procedure yields excellent recovery 
MSE. However, for "heavy-tailed" signals, such as those gen- 
erated from Student's t-distributions, our experience indicates 
that the EM-learned values of {9k} tend to gravitate towards 
the outliers in {x n }^ =1 , resulting in an overriding of px (•) and 
thus poor reconstruction MSE. For such heavy-tailed signals, 



(58) 



We start by setting ^-J(u>,£) = 0, which yields 



N 



/ Px(x n ;q )Af(x n ;r n ,n n ) d , / . s e 



(59) 



N 

- E 

71=1 



Px (x n ; q l )Af{xn \r„,n r n ) \W{x n ;6' k ,(t>l) 
C(r n ,^ , n ;q i ) Px(x n ;u>,q\ u ) 



(60) 
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Like in d48l and (l54l ), the above integral is difficult to evaluate, 
and so we approximate u> « w*, which reduces the previous 
equation to 



N 

£ 

ra=l 



(61) 



Multiplying both sides by u> l k for k = 1 , . . . , L, summing over 
k, employing the fact 1 = Efc w fc' anc ^ simplifying, we obtain 
the equivalent condition 



N 



n=l 
AT 



^Pr{x„^0|y;q 1 }. 



(62) 



(63) 



Plugging d63l into ( |6TT > and multiplying both sides by Wfe, the 
derivative-zeroing value of Uk is seen to be 

Y,n=iL n xicj kN(x n ; Oi,<t>i)JV(x n ;r n , ^ r n )/((f n , fa q l ) 



EliPr{x„^0|y;^} 



(64) 



where, if we use u>k ~ on the right of d64l i. then we obtain 



En=l Pr { x « 7^ 0, fc„ = fc | y; g 4 } 



EliPrR^Oly;^} 

Note that the numerator and denominator of d65l l can be 
computed from GM-GAMP outputs via dBll ), ( fT6l ), and Q. 

Although, for the case of GM priors, approximations were 
used in the derivation of the EM updates d50l l. d56i l. and d65l ). 
it is interesting to note that, in the case of L = 1 mixture 
components, these approximate EM-GM updates coincide 
with the exact EM-BG updates derived in Section IIII-BI In 
particular, the approximate-EM update of the GM parameter 
9\ in (f50T > coincides with the exact-EM update of the BG 
parameter 9 in (l35l l. the approximate-EM update of the GM 
parameter <fii in (l56l l coincides with the exact-EM update of 
the BG parameter in (HTt . and the approximate-EM update 
of the GM parameter lu\ in (l(55l l reduces to the fixed value 1. 
Thus, one can safely use the GM updates above in the BG 
setting without any loss of optimality. 

D. EM Initialization 

Since the EM algorithm is guaranteed to converge to only a 
local maximum of the likelihood function, proper initialization 
of the unknown parameters q is essential. Here, we propose 
initialization strategies for both the "sparse" and "heavy- 
tailed" modes of operation, for a given value of L. 

For the "sparse" mode, we set the initial sparsity rate 
A equal to the theoretical noiseless Lasso PTC, i.e., A° = 
M pPse(# ), where EU) 



(65) 



N ' 



PSe(^) = max c>0 



l-2£[(l + c 2 )$( C )- C( /,(c)] 



(66) 



l-c 2 -2[(l + c 2 )$(c)-c0(c)] 

describes the maximum value of jj supported by Lasso for a 
given j^r, and where $(•) and <f>(-) denote the cdf and pdf of 
the Af(0, 1) distribution, respectively. Using the energies \\y\\\ 



and and an assumed value of SNR°, we initialize the 

noise and signal variances, respectively, as 

Mif)° 



lylll 



(SNR° + 1)M 



Ivlli 



|A||f,A° 



(67) 



where, in the absence of (user provided) knowledge about 
the true SNR = || Ax|||/|H||, we suggest SNR° = 100. 
Then, we uniformly space the initial GM means 9° over 
\ ~2L 1 i ^Xr]> an< ^ subsequently fit the mixture weights u>° 
and variances 0° to the uniform pdf supported on [—0.5, 0.5] 
(which can be done offline using the standard approach to 
EM-fitting of GM parameters, e.g., 12T1 p. 435]). Finally, we 
multiply 0° by y / !2ip° and 0° by 12ip° to ensure that the 
resulting signal variance equals ip°. 

For the "heavy-tailed" mode, we initialize A and ip° as 
above and set, for fe = 1, . . . , L, 



"1 



J^{\\y\\l-M^) 
VI" 



\Af F \* 



and 9l = 0. (68) 



E. Selection of GM Model Order L 

Until now, the GM model order L has been treated as fixed 
and known. Indeed, one approach to model-order selection is 
to choose a fixed value of L that is thought to be large enough 
to model the essential structure within the true signal pdf and, 
after an appropriate initialization of the 3L + 2 parameters 
in q (such as that described in Section IHI-Dl i. to apply the 
EM-GM- AMP algorithm summarized in Table [II] 

An alternative approach to model-order selection is to start 
with the model order L = 1 (i.e., Bernoulli-Gaussian px{-)) 
and increment L in steps of one, stopping as soon as negligible 
benefits are observed (e.g., \\xl — ^i-illi/ll^i-ilii < tol) 
or a predefined L max has been reached. Here, EM-GM-AMP 
would be re-run as described above for each new value of L. 
This "greedy" approach would relieve the user of the task of 
choosing an appropriate L and initializing the corresponding 
3L + 2 parameters in q. In the remainder of this section, we 
propose a particular instantiation of this greedy approach. 

When growing the model-order from L to L + 1, we 
choose to split the mixture component fc* E {1, . . . , L} with 
the "worst fit" into two new components. As a criterion 
for the worst-fitting mixture component, one could use local 
Kullback-Leibler divergence, as in Ueda's split-and-merge-EM 
algorithm f23l , or similar. Using fc* to denote the index of 
the mixture-component-to-split, the subset of signal coefficient 
indices n that are a-posteriori most-probably associated with 
k* is identified, i.e., 

9tfc„ — {n^UX: argmaxPr{x„^0, k n = k \ y; qr} = fc* } ,(69) 

k 

where 9T = {n : Pr{x n ^ | y; q} > 0.5} is the subset of 
coefficient indices that are a-posteriori most-probably non-zero 
(i.e., the GAMP support estimate) and where q denotes the 
current estimate of the parameters. To simplify the notation, 
we henceforth assume, without loss of generality, that fc* = L. 

To split the L th mixture component, we then replace its 
mean 6l, variance <j>^, and weight cj^ with two new values for 
each (i.e., 6*" ew and would replace 8^), resulting in a new 
parameter vector q new containing 3(L+l)+2 terms. Moreover, 
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L+l 



y L + l — ' 

= uj l /2; 

and 
, .new 



6 new 



new 

L + l 



rather than considering only a single possibility for q new , we 
allow that S > 1 hypotheses {q" ew }f =1 are considered, where 
the values in each q" ew are produced by some variation on the 
following two strategies: 

1) Split-mean(a): 6f w = 6 L - a, =6 L +a, 

= <t> L , and w2 ew = w 

2) Split-variance(6): <p£» 
e new = gnaw = ^ and ^new = w new = Wj ./ 2> 

where a, 6 > are design parameters. Note that, by consid- 
ering several distinct values of a and/or b, we get S > 2. 
Finally, to choose among the hypothesized splits {<7s ew }f =1 , 
one could, e.g., select the one that maximizes the likelihood 
E{lnp x ( X ;q"™)\y;q}. 

We investigated the performance of EM-GM-AMP under 
this "greedy" model-order selection procedure in 0. Com- 
paring the numerical results in |3] (where greedy model- 
order selection was used for the "sparse" mode) to those in 
Section |IV] of this work (where a fixed model-order of L = 3 
was employed for the sparse mode), it is evident that the 
greedy scheme from (3j performs slightly better. However, that 
greedy scheme has a significantly higher complexity, mainly 
because the EM-GM-AMP algorithm in Table [II] must be re- 
run for each tested value of L. For this reason, we focus on 
the fixed model-order approach in the remainder of this work, 
and regard accurate and fast model-order-selection as a topic 
of future study. 

F. EM-GM-AMP Summary and Demonstration 

The fixed-L EM-GM-AMlH algorithm developed in the 
previous sections is summarized in Table IHl For EM-BG-AMP 
(as previously described in iff)), one would simply run EM- 
GM-AMP with L = 1. 

To demonstrate EM-GM-AMP's ability to learn the un- 
derlying signal distribution, Fig. [T] shows examples of the 
GM-modeled signal distributions learned by EM-GM-AMP 
in both "sparse" and "heavy-tailed" modes. To create the 
figure, we first constructed the true signal vector x £ M. N 
using N = 2000 independent draws of the true distribution 
Px(') shown in each of the subplots. Then, we constructed 
the matrix A £ M. MxN and the noise vector w £ M. AI using 
independent zero-mean Gaussian draws, with M = 1000 and 
the noise variance adjusted to achieve SNR=25 dB in the 
resulting observations y = Ax + w. Finally, we ran EM- 
GM-AMP according to Table HH and plotted the GM approx- 
imation px{x;q l ) from using the learned pdf parameters 
q l = [X\uj\e\4>\ Figure □ confirms that EM-GM- 
AMP is successful in learning a reasonable approximation 
of the unknown true pdf px(-) from the noisy compressed 
observations y, in both sparse and heavy-tailed modes. 

IV. Numerical Results 

In this section we report the results of a detailed numerical 
study that investigate the performance of EM-GM-AMP under 

5 Matlab code available at http://www.ece.osu.edu/~schniter/EMturboGAMP 
Empirically, we have observed that the EM update for tp works better 
with the /if n term in (23) weighted by and suppressed until later EM 
iterations. We conjecture that this is due to bias in the finite-sample GAMP 
variance estimates /jtfL. 
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Fig. 1. True and EM-GM-AMP-learned versions of the signal distribution 
Px{ x ) = ^fx( x ) + (1 — X)S(x). The top subplot shows "sparse" mode 
EM-GM-AMP run using GM-order L = 3 on a sparse signal whose non- 
zero components were generated according to a triangular mixture, whereas 
the bottom subplot shows "heavy-tailed" EM-GM-AMP run using L = 4 on a 
Student's-t signal with rate parameter q = 1.67 (defined in i70\ ). The density 
of the continuous component Xfx ( x ) is marked on the left axis, while the 
mass of the discrete component (1 — \)S(x) is marked on the right axis. 



Initialize L as described in Section IIII-EI 




Initialize q° as described in Section IIII-DI 




Initialize x = 0. 




for i = 1 to 7 max do 




Generate x\ z\ (fJ, z )\ tt\ {^l^k^Uk-l 


using GM-GAMP 


with (f -1 (see Tabled. 




if \\x % — a; I_1 ||2 < r e m||i l_1 1|| then 




break. 




end if 




Compute A 1 from 7r I_1 as described in )29t. 




for k = 1 to L do 




if sparse mode enabled then 




Compute 65 from 7r l 1 , y\ , {/3! as described in 


(50). 




else if heavy-tailed mode enabled then 




Set 0* = 0. 




end if 




Compute 4>i from tt*" 1 , 7^ 1 , i/*" 




described in (56). 


Compute ijj 1 from 7r l_1 and {f3 1 l ~ 1 }f'_ 1 as 


described in (65). 


end for 




Compute t/i 1 from z l and (fJ, z Y as in (23)0 




end for 





TABLE II 

THE EM-GM-AMP ALGORITHM (FIXED-L CASE) 

both noiseless and noisy settings. For all experiments, we set 
the GM-GAMP tolerance to Tg am p = 10 -5 and the maximum 
GAMP-iterations to T max = 20 (recall Table HJ), and we set the 
EM tolerance to T em = 10~ 5 and the maximum EM-iterations 
to /max = 20 (recall Table Hill. In "sparse" mode, we use GM 
order L = 3, while in "heavy-tailed" mode, we use L = 4. 

A. Noiseless Phase Transitions 

We first describe the results of experiments that computed 
noiseless empirical phase transition curves (PTCs) under 
three sparse-signal distributions. To evaluate each empirical 
PTC, we fixed N = 1000 and constructed a 30 x 30 grid 
where (M, K) were chosen to yield a uniform sampling 
of oversampling ratios M- £ [0.05, 0.95] and sparsity ratios 
jfe £ [0.05,0.95]. At each grid point, we generated R = 100 
independent realizations of a A'-sparse signal x from a spec- 
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ified distribution and an M x N measurement matrix A with 
i.i.d Af(0, M _1 ) entries. From the noiseless measurements 
y = Ax, we recovered the signal x using several algorithms. 
A recovery x from realization r £ {1, ... ,R} was defined 
a success if the NMSE = \\x - x\\%/\\x\\l < 10" 4 , and the 
average success rate was defined as S = Y^=i where 
S r = 1 for a success and S r = otherwise. The empirical 
PTC was then plotted, using Matlab's contour command, 
as the S = 0.5 contour over the sparsity-undersampling grid. 

Figures [2]|4] show the empirical PTCs for four recovery 
algorithms: the proposed EM-GM-AMP algorithm (in "sparse" 
mode), the proposed EM-BG-AMP algorithm, a genie-tunecQ 
GM-AMP that uses the true parameters q = [A, u>, 0, </>, i/'], 
and the Donoho/Maleki/Montanari (DMM) Lasso-style AMP 
from ifTUI . For comparison, Figs. [2]|4] also display the theoreti- 
cal Lasso PTC ( l66l ). The signals were generated as Bernoulli- 
Gaussian (BG) in Fig. [2] (using mean 8 — and variance 
0=1 for the Gaussian component), as Bernoulli in Fig. [3] 
(i.e., all non-zero coefficients set equal to 1), and as Bernoulli- 
Rademacher (BR) in Fig. [4] (i.e., non-zero coefficients chosen 
uniformly at random from the set { — 1,1}). 

For all three signal types, Figs. [2]|4] show that the empirical 
PTC of EM-GM-AMP significantly improves on the empirical 
PTC of DMM-AMP as well as the theoretical PTC of Lasso. 
(The latter two are known to converge in the large system 
limit Q0|.) For BG signals, Fig. shows that EM-GM-AMP 
and EM-BG-AMP both yield PTCs that are nearly identical 
to that of genie-GM-AMP, suggesting that our EM-learning 
procedure is working well. Note that, for these BG signals, 
the model employed by EM-BG-AMP is perfectly matched to 
the true signal, whereas that employed by EM-GM-AMP (with 
L = 3) has the potential to over-fit the true signal, although 
Fig. |2 suggests that this does not happen. For Bernoulli signals, 
Fig. [3] shows EM-BG-AMP performing nearly the same as 
genie-GM-AMP, and EM-GM-AMP performing even better 
than genie-GM-AMP. The latter behavior is due to EM-GM- 
AMP's ability to do per-realization parameter tuning, whereas 
genie-GM-AMP employs the best set of fixed parameters over 
all realizations. Finally, for BR signals, Fig. [4] shows EM-GM- 
AMP performing significantly better than EM-BG-AMP, since 
the former has the ability to accurately model the BR distri- 
bution (with L > 2 mixture components), whereas the latter 
(with a single mixture component) does not. Figure 0] again 
shows EM-GM-AMP performing slightly better than genie- 
GM-AMP in some regions of the sparsity-undersampling plane 
as a result of per-realization parameter tuning. 

B. Noisy Sparse Signal Recovery 

Figures |5]|7] show NMSE for noisy recovery of BG, 
Bernoulli, and BR signals, respectively. To construct these 
plots, we fixed N = 1000, K = 100, SNR = 25 dB, 
and varied M. Each data point represents NMSE averaged 
over R = 500 realizations. For comparison, we show the 
performance of the proposed EM-GM-AMP (in "sparse" 

7 For genie-tuned GM-AMP, for numerical reasons, we set the noise variance 
at ij> = 10 -6 and, with Bernoulli and BR signals, the mixture variances at 

<t>k = io- 2 . 
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Fig. 2. Empirical PTCs and Lasso theoretical PTC for noiseless recovery of 
Bernoulli-Gaussian signals. 




0.4 0.5 0.6 

M/N 

Fig. 3. Empirical PTCs and Lasso theoretical PTC for noiseless recovery of 
Bernoulli signals. 




theoretical Lasso 
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Fig. 4. Empirical PTCs and Lasso theoretical PTC for noiseless recovery of 
Bernoulli-Rademacher signals. 



mode), EM-BG-AMP, genie-tuned^ Orthogonal Matching Pur- 
suit (OMP) J24), genie-tunecP Subspace Pursuit (SP) 1251 , 
Bayesian Compressive Sensing (BCS) |fl9l . Sparse Bayesian 
Learning lfT8l (via the more robust T-MSBL ll26l ). de-biased 



8 We ran both OMP (using the implementation from 
http://sparselab.stanford.edu/OptimalTuning/code.htm I and SP under 10 
different sparsity assumptions, spaced uniformly from 1 to 2K, and reported 
the lowest NMSE among the results. 
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genie-tunecfl Lasso (via SPGL1 112710 , and Smoothed-£ (SLO) 
11281 . All algorithms were run under the suggested defaults, 
with , noise=small' in T-MSBL. 

For BG signals, Fig. shows that EM-GM-AMP and EM- 
BG-AMP together exhibit the best performance among the 
tested algorithms, reducing the M/N breakpoint (i.e., the lo- 
cation of the knee in the NMSE curve, which represents a sort 
of phase transition) from 0.3 down to 0.26, but also improving 
NMSE by ps 1 dB relative to the next best algorithm, which 
was BCS. For Bernoulli signals, Fig. [6] shows much more 
significant gains for EM-GM-AMP and EM-BG-AMP over the 
other algorithms: the M/N breakpoint was reduced from 0.4 
down to 0.32, and the NMSE was reduced by « 8 dB relative 
to the next best algorithm, which was T-MSBL in this case. 
Finally, for BR signals, Fig. [7] shows a distinct advantage for 
EM-GM-AMP over the other algorithms, including EM-BG- 
AMP, due to the former's unique ability to accurately model 
the BR signal prior. In particular, for M /N > 0.36, EM-GM- 
AMP reduces the NMSE by at least 8 dB relative to the best 
of the other algorithms (which was either EM-BG-AMP or 
T-MSBL depending on the value of M/N) and reduces the 
M/N breakpoint from 0.38 down to 0.36. 
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Fig. 5. NMSE versus undersampling ratio M/N for noisy recovery of 
Bernoulli-Gaussian signals. 

To investigate each algorithm's robustness to AWGN, we 
plotted the NMSE attained in the recovery of BR signals 
with N = 1000, M = 500, and K = 100 as a function 
of SNR in Fig. [8] where each point represents an average 
over R = 100 problem realizations. All algorithms were under 
the same conditions as those reported previously, except that 
T-MSBL used , noise=small' when SNR > 22dB and 
, noise=mild' when SNR < 22 dB, as recommended in 
11291 . From Fig. [S] we see that the essential behavior observed 
in the fixed-SNR BR plot Fig. [7] holds over a wide range 
of SNRs. In particular, Fig. shows that EM-GM-AMP 
yields significantly lower NMSE than all other algorithms over 
the full SNR range, while EM-BG-AMP and T-MSBL yield 
the second lowest NMSE (also matched by BCS for SNRs 
between 30 and 40 dB). Note, however, than T-MSBL must be 
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Fig. 6. NMSE versus undersampling ratio M/N for noisy recovery of 
Bernoulli signals. 
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Fig. 7. NMSE versus undersampling ratio M/N for noisy recovery of 
Bernoulli-Rademacher signals. 



given some knowledge about the true noise variance in order 
to perform well |29), unlike EM-GM-AMP or EM-BG-AMP. 
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hypothesized tolerances a' 2 S {0.1,0.2,.. 
lowest NMSE among the results. 



s.t. \\y — Ax || 2 < <J, for 
, 1.5} X Mip, and reported the 
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Fig. 8. NMSE versus SNR for noisy recovery of Bernoulli-Rademacher 
signals. 



C. Heavy-Tailed Signal Recovery 

In many applications of compressive sensing, the signal to 
be recovered is not perfectly sparse, but rather "heavy tailed," 
in that there are a few large coefficients and many small (but 
nonzero) ones. To investigate algorithm performance in this 
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setting, we constructed the signal vector x as i.i.d Student's-t, 
i.e., with prior pdf 



Px(x;q) 



a r((g+i)/2)) 

v / 2¥r(g/2) 



-(9+l)/2 



(70) 



under the (non-compressible) rate q = 1.67, which has been 
shown to be an excellent model for wavelet coefficients of 
natural images l30l . For such signals, Fig. [9] plots NMSE 
versus the number of measurements M for fixed N = 1000, 
SNR = 25 dB, and an average of R = 500 realizations. 
Figure [9] shows EM-GM-AMP (here run in "heavy-tailed" 
mode) outperforming all other algorithms under test0 We 
have also verified (in experiments not shown here) that "heavy- 
tailed" EM-GM-AMP exhibits similarly good performance 
with other values of the Student's-t rate parameter q. 

It may be interesting to notice that, with the perfectly 
sparse signals examined in Figs. |5]|7] the mixed-norm ap- 
proaches (i.e., SL0 and SPGL1) performed relatively poorly, 
the relevance-vector-machine (RVM)-based approaches (i.e., 
BCS, T-MSBL) performed relatively well, and the greedy 
approaches (OMP and SP) performed in-between. Meanwhile, 
with the heavy-tailed signals in Fig. [9] the situation was re- 
versed: the mixed-norm approaches performed relatively well, 
whereas the RVM approaches performed relatively poorly. 
Thus, it appears that the mixed-norm approaches are somehow 
better "tuned" to heavy-tailed signals, whereas the RVM 
approaches are better tuned to sparse signals. For all signal 
types, though, the best recovery performance came from EM- 
GM-AMP, and we attribute this behavior to EM-GM-AMP's 
ability to tune itself to the signal realization at hand. 
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Fig. 9. NMSE versus undersampling ratio M/N for noisy recovery of 
Student's-t signals. 



D. Runtime and Complexity Scaling with N 

Next we investigated how complexity scales with signal 
length N by evaluating the runtime of each algorithm on a 
typical personal computer. For this, we fixed K/N = 0.1, 
M/N = 0.5, SNR = 25 dB and varied the signal length 
N. Figure [10] shows the runtimes for noisy recovery of a 



Bernoulli-Rademacher signal, while Fig. QT| shows the cor- 
responding NMSEs. In these plots, each datapoint represents 
an average over R = 50 realizations. The algorithms that we 
tested are the same ones that we described earlier. However, to 
fairly evaluate runtime, we configured some a bit differently 
than before. In particular, for genie-tuned SPGL1, in order 
to yield a better runtime-vs-NMSE tradeoff, we reduced the 
tolerance grid (recall footnote|9]l to a 2 <E {0.6, 0.8, . . . , 1.4} x 
Mip and turned off debiasing. For OMP and SP, we used 
under the fixed support size i^iasso = Mpse(j^) rather than 
searching for the size that minimizes NMSE over a grid 
of 10 hypotheses, as before. Otherwise, all algorithms were 
run under the suggested defaults, with T-MSBL run under 
, noise=small' and EM-GM-AMP run in "sparse" mode. 

The complexities of EM-BG-AMP and EM-GM-AMP are 
dominated by one matrix multiplication by A and A J per 
iteration, and the number of iterations is invariant to M and 
N, Thus, when these matrix multiplications are explicitly 
implemented and A is dense, the total complexity of EM- 
BG-AMP and EM-GM-AMP should scale as 0{MN). This 
scaling is indeed visible in the runtime curves of Fig. [Tol 
There, 0(MN) becomes 0{N 2 ) since the ratio M/N was 
fixed, and the horizontal axis plots N on a logarithmic scale, 
so that this complexity scaling manifests, at sufficiently large 
values of N, as a line with slope 2. Figure [10] confirms that 
genie-tuned SPGL1 also has the same complexity scaling, 
albeit with longer overall runtimes. Meanwhile, Fig. [Tol shows 
T-MSBL, BCS, SL0, OMP, and SP exhibiting a complexity 
scaling of 0(N 3 ) (under fixed K/N and M/N), which results 
in orders-of-magnitude larger runtimes for long signals (e.g., 
N > 10 4 ). With short signals (e.g., N < 1300), though, 
OMP, SP, SL0, and SPGL1 are faster than EM-GM-AMP. 
Finally, Fig. QT| verifies that, for most of the algorithms, the 
NMSEs are relatively insensitive to signal length N when the 
undersampling ratio M/N and sparsity ratio K/M are both 
fixed, although the performance of EM-GM-AMP improves 
with N (which is not surprising in light of AMP's large- 
system-limit optimality properties ITT31 ) and the performance 
of BCS degrades with N. 

The proposed EM-GM-AMP and EM-BG-AMP algorithms, 
as well as SPGL1, can also handle the case where multipli- 
cation by A and A J is implemented using a fast algorithm 
like the fast Fourier transform (FFT0 which reduces the 
complexity to 0(N log N), and avoids the need to store A in 
memory — a serious problem when MN is large. The dashed 
lines in Figs. [TolfTTI (labeled "f ft") show the average runtime 
and NMSE of EM-GM-AMP, EM-BG-AMP, and SPGL1 in 
case that A was a randomly row-sampled FFT As expected, 
the runtimes are dramatically reduced. While EM-BG-AMP 
retains its place as the fastest algorithm, SPGL1 now runs 
1.5 x faster than EM-GM-AMP (at the cost of 14 dB higher 
NMSE). 



In this experiment, we ran both OMP and SP under 10 different sparsity 
hypotheses, spaced uniformly from 1 to A"| ass0 = Mpqe(M-), and reported 
the lowest NMSE among the results. 



1 1 For our FFT-based experiments, we used the complex- valued versions of 
EM-BG-AMP, EM-GM-AMP, and SPGL1. 
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Fig. 10. Runtime versus signal length N for noisy recovery of Bernoulli- 
Rademacher signals. 
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Fig. 11. NMSE versus signal length N for noisy recovery of Bernoulli- 
Rademacher signals. 



E. Example: Compressive Recovery of Audio 

As a practical example, we experimented with the recovery 
of an audio signal from compressed measurements. The full 
length-81920 audio signal was first partitioned into T blocks 
{u t }f =1 of length N. Noiseless compressed measurements 
y t = <&u t € M M were then collected using M = N/2 
samples per block. Rather than reconstructing u t directly 
from y t , we first reconstructec(3 the transform coefficients 
x t = * T Mt, using the (orthogonal) discrete cosine transform 
(DCT) * G M ArxJV , and later reconstructed u t via u t = ^x t . 
Our effective sparse-signal model can thus be written as 
y t = Ax t with A = 4"$'. We experimented with two types of 
measurement matrix <&: i.i.d Gaussian and random selection 
(i.e., containing rows of the identity matrix selected uniformly 
at random), noting that the latter allows a fast implementation 
of A and A J . Table [HI] shows the resulting time-averaged 
NMSE, i.e., TNMSE 4 i £f=i ||it t -M t || 2 /||M t || 2 , and total 
runtime achieved by the previously described algorithms at 
block lengths N = 1024, 2048, 4096, 8192, which correspond 
to T = 80,40,20,10 blocks, respectively. The numbers 

12 Although one could exploit additional structure among the multiple- 
timestep coefficients {&t}T=i for improved recovery (e.g., sparsity clustering 
in the time and/or frequency dimensions, as well as amplitude correlation in 
those dimensions) as demonstrated in 1311 , such techniques are outside the 
scope of this paper. 



reported in the table represent an average over 50 realizations 
of <I>. For these experiments, we configured the algorithms 
as described in Section IIV-CI for the heavy-tailed experiment 
except that, for genie-SPGLl, rather than using ip = 0, we 
used tp = 10 -6 for the tolerance grid (recall footnote [9]) 
because we found that this value minimized TNMSE and, for 
T-MSBL, we used the setting prune_gamma = 10 -12 as rec- 
ommended in a personal correspondence with the author. For 
certain combinations of algorithm and blocklength, excessive 
runtimes prevented us from carrying out the experiment, and 
thus no result appears in the table. 

Table [Til] shows that, for this audio experiment, EM-GM- 
AMP and SL0 performed best in terms of TNMSE, and EM- 
BG-AMP performed second best. As in the synthetic exam- 
ples presented earlier, we attribute EM-GM-AMP's excellent 
TNMSE to its ability to tune itself to whatever signal is 
at hand. As for SLO's excellent TNMSE, we reason that 
it had the good fortune of being particularly well-tuned to 
this audio signal, given that it performed relatively poorly 
with the signal types used for Figs. |5]|8] From the runtimes 
reported in Table UTTI we see that, with i.i.d Gaussian <& and the 
shortest block length (N = 1024), OMP is by far the fastest, 
whereas EM-BG-AMP and EM-GM-AMP are the slowest. 
But, as the block length grows, EM-BG-AMP and EM-GM- 
AMP achieve better and better runtimes as a consequence 
of their excellent complexity scaling, and eventually become 
the fastest of the algorithms under test (as shown with i.i.d 
Gaussian 3? at N = 8192). For this audio example, the large- 
block regime may be the more important, because that is 
where all algorithms give smallest TNMSE. Next, looking 
at the runtimes under random-selection <fr, we see dramatic 
speed improvements for EM-GM-AMP, EM-BG-AMP, and 
SPGL1, which were all able to leverage Matlab's fast DCT. 
In fact, the total runtimes of these three algorithms decrease 
as iV is increased from 1024 to 8192. We conclude by noting 
that EM-BG-AMP (at N = 8192 with random selection *) 
achieves the fastest runtime in the entire table while yielding 
a TNMSE that is within 1 dB of the best value in the entire 
table. Meanwhile, EM-GM-AMP (at N = 8192 with random 
selection $) yields the best TNMSE in the entire table while 
taking only about twice as long to run as the fastest time in 
the entire table. 

V. Conclusions 

Those interested in practical compressive sensing face the 
daunting task of choosing among literally hundreds of signal 
reconstruction algorithms (see, e.g., [32]). In testing these 
algorithms, they are likely to find that some work very well 
with particular signal classes, but not with others. They are 
also likely to get frustrated by those algorithms that require the 
tuning of many parameters. Finally, they are likely to find that 
some of the algorithms that are commonly regarded as "very 
fast" are actually very slow in high-dimensional problems. 
Meanwhile, those familiar with the theory of compressive 
sensing know that the workhorse Lasso is nearly minimax 
optimal, and that its phase transition curve is robust to the 
nonzero-coefficient distribution of sparse signals. However, 
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COMPRESSIVE AUDIO RECOVERY. 



they also know that, for most signal classes, there is a large 
gap between the MSE performance of Lasso and that of the 
MMSE estimator derived under full knowledge of the signal 
and noise statistics fTTl . Thus, they may wonder whether there 
is a way to close this gap by designing a signal reconstruction 
algorithm that both learns and exploits the signal and noise 
statistics. 

With these considerations in mind, we proposed an em- 
pirical Bayesian approach to compressive signal recovery 
that merges two powerful inference frameworks: expectation 
maximization (EM) and approximate message passing (AMP). 
We then demonstrated — through a detailed numerical study — 
that our approach, when used with a flexible Gaussian-mixture 
signal prior, achieves a state-of-the-art combination of recon- 
struction error and runtime on a very wide range of signal 
types in the high-dimensional regime. 
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